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In this talk we show how the sign problem, occurring in dynamical simulations of random matri- 
ces at nonzero chemical potential, can be avoided by judiciously combining matrices into subsets. 
One can prove that these subsets have real and positive weights such that importance sampling 
can be used in Monte Carlo simulations. The number of matrices per subset is proportional to the 
(~| matrix dimension. We measure the chiral condensate and observe that the statistical error is inde- 

pendent of the chemical potential and grows linearly with the matrix dimension, which contrasts 
strongly with its exponential growth in reweighting methods. 
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1. Introduction 

Dynamical Monte Carlo simulations of QCD are seriously hampered at nonzero chemical 
potential /i because the fermion determinant becomes complex, causing the notorious sign problem 
[1]. The sign problem in QCD can be explored using random matrix theory (RMT) [2] because 
of the equivalence between QCD in the £-regime and RMT [3]. Although many observables in 
unquenched RMT have been computed analytically, it is interesting to investigate if one could also 
access these numerically. Dynamical simulations of random matrices at nonzero jU also suffer from 
a sign problem and can therefore be used as a playground for algorithmic developments. In this talk 
we will present a subset method, which solves the sign problem in the RMT case by judiciously 
combining matrices into subsets with real and positive fermionic weights [4]. 



2. Fermion determinant and sign problem in QCD 

After integration over the fermion fields the QCD partition function can be written as: 

Nf 



Zqcd= / e-^'^Y['^Qt[D^,^^] , (2.1) 

/=i 

^ V ' 

MCMC weight function P ? 

where only the integration over the gauge fields remains, Sq is the gauge action and is 
the Dirac operator for a quark of mass nif at chemical potential jj.. As long as the weight factors 
are real and positive, the functional integral in lattice QCD is evaluated by a Markov chain Monte- 
Carlo (MCMC) simulation using importance sampling, and the expectation value of an observable 
O is approximated by the sample average of the measurements on A^mc configurations: 

The fermion determinant det[D^ „,] is real and positive for jJ. = 0, but becomes complex when 
jJ- ^0. In the latter case the fermion determinant can no longer be interpreted as a probabilistic 
weight in MCMC simulations and we are confronted with the sign problem. Methods to perform 
measurements at finite chemical potential, by circumventing the sign problem, generically require 
a computing time which grows exponentially with the volume. This is, for example, the case in 
reweighting methods, where the ensemble is sampled according to an auxiliary weight function 
and the results aie reweighted appropriately. 



3. Random matrix theory 

In the £-regime, QCD is equivalent to chiral random matrix theory, both at zero and nonzero 
chemical potential [3]. In the two-matrix model of Osborn [5] the random matrices 0i and 02 are 
complex {N + v) xN matrices distributed according to the unquenched partition function 

/• 

Zy^{H;{mf})= d(l)id(l)2wi^i)w{^)Y[<i^iD^^^^{^i,^) (3.1) 

/=i 
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sign problem 



Figure 1: Average phase factor (e^'^) as a function of the chemical potential a = IN^l^ for m = (left) and 
as a function of the mass m = 2Nm for a = 1 (right). A small average phase factor corresponds to a strongly 
fluctuating phase and is evidence of the sign problem. 

with Gaussian weights w(0) = {N /n)'^^'^^^^ exp{—Ntr^'' ^) and Nf dynamical quarks of masses 
my at a chemical potential /i, whose Dirac operator (with v zero modes) is given by: 

The dynamics of the random matrix model crucially depend on the determinant of the Dirac op- 
erator. Just as in QCD, the Dirac matrix D^ ,„ is non-Hermitian in this random matrix model for 
/X 7^ 0: Its determinant is complex and can be written as det[D^ ,„] = Re^^ . The average phase factor 
(e^'^) reflects the fluctuations of the fermion determinant and characterizes the strength of the sign 
problem in dynamical simulations. It was computed analytically in refs. [2, 6, 7] and is shown in 
fig. 1, where we highlighted the parameter region where the sign problem occurs. 

In this paper we examine the sign problem in dynamical simulations of chiral random matrices. 
Before presenting a new solution for the sign problem we briefly describe the standard reweighting 
method, which we will later use to compare with the new results and to verify the onset of the sign 
problem. 



4. Reweighting 

The reweighting method can be used to circumvent the sign problem and perform measure- 
ments in dynamical simulations with complex weights. The ensemble average of an observable 
y{x) in an ensemble with weight w{x) is defined by 

J dx w[x) 

In the reweighting method one introduces an auxiliary ensemble with weight Waux(-^) and 
rewrites the previous equation as 

(^)_,^J''^-'---^W^^W,(^^L. (4.2, 
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If the ensemble Waux is chosen to be real and positive it can be sampled using importance sampling 
methods and the result of eq. (4.2) can be evaluated in a Monte Carlo simulation. Typical examples 
for Waux are the quenched, phase-quenched, -quenched and sign-quenched ensembles. 

The problem with reweighting methods is that the work needed to make reliable measurements 
on the statistical ensemble grows exponentially with volume and chemical potential because it 
involves the computation of exponentially small reweighting factors from a statistical sampling of 
largely canceling contributions [1]. 



5. Subset method for dynamical RMT simulations 

Below we describe a subset method which solves the sign problem for dynamical simulations 
of the Osbom model and was first introduced in ref. [4]. 

For any given random matrix pair = (0i , ^2) we introduce a set of matrices 

= I e,):dn = —An = 0,... ,A^,-1 J , (5.1) 

containing A'^, orthogonal rotations i/a = (i/^i , i//2) of defined by 



Yi{^;6)\ _ I COS0 sin0 W (/>i 
Y2{(l>;0) I -sine COS0/ I 



(5.2) 



The subset construction (5.1) allows for a partial resummation of the original random matrix parti- 
tion function (3.1), which can be rewritten as an equivalent partition function over subsets ^l: 



j dQ.W{Q.)oci{^l,m). (5.3) 



The subset weights were factorized in a Gaussian part W{Q.) = w{\\f\ {(j);d))w{Y2{<P',Q)), which is 
independent of 6 because of the orthogonal rotations (5.2), and a fermionic weight 



Ga{n,m)= £ det^/D^,„,(v^((/.;0„)), (5.4) 

n=0 



which is a sum of complex determinants. The equivalence of the partition functions (3.1) and 
(5.3) rests on the observation that there is a subset Q.{<j)) for each configuration (j) = ^2) of the 
original partition function, so that the set of all subsets forms an A(s-fold covering of the original 
RMT ensemble. 

The subset method solves the sign problem because of the following positivity theorem: For 
any subset Q. given by eq. (5.1) the fermionic subset weight OQ,{pL,m) is real and positive if 
Ns > NfN (for arbitrary m and /i < 1 ). 

This theorem results from the following identity, which relates the fermionic subset weights at 
nonzero and zero chemical potential. For arbitrary jj. and m: 

aU^,m) = {l-^^rf(^+'i)aa U y=^] > (5-5) 
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for any Q. constructed according to (5.1) /f > NfN. The proof of this identity will be given in a 
forthcoming publication. From this identity the positivity of an(/i,m) is easily derived: For /i = 
the determinants of all the matrices in the subset are real and positive, as all the eigenvalues of 
the Dirac matrix come in complex conjugate pairs in this case (for arbitrary real mass). Therefore, 
eq. (5.5) implies that the fermionic weight OQ,{pL,m) is real and positive for /i < 1. Moreover, for 
jLt = 1 and m = the sum of determinants is exactly zero, which corresponds to the case of maximal 
non-hermiticity. 

Note that eq. (5.5) is an extension to arbitrary mass of the identity originally given in ref. [4], 
which only covered the massless case, while an inequality described the case m ^ 0. 

6. Simulations 

The positive subset weights W{D.)aQ,{lJ.,m) were used to generate subsets of random matrices 
and sample the partition function (5.3) with a Metropolis algorithm. In practice the subset size is 
set to Ns = NfN + 1, which is the minimum value for which the positivity of the fermionic weights 
is guaranteed. Successive subsets in the Markov chain are generated as follows: 

• randomly choose a configuration in the current subset, 

• generate a new configuration by making a random step, 

• construct the subset corresponding to the new configuration, 

• apply the accept-reject step to the proposed subset using the positive subset weights. 

This algorithm satisfies detailed balance and ergodicity such that the partition function will be sam- 
pled correctly by the MCMC algorithm. In the subset method the sample average O measured over 
a sample of A^mc subsets H^, approximating the ensemble average in the original RMT ensemble, 
is computed by 

O.^nr = I I ^ ^ 0,.Uy'% (6.1) 

where G Q.^ and one takes into account that the matrices inside the subsets yield different values 
for the measured observable. 

7. Results 

We applied the subset method to compute the chiral condensate 

in the RMT model. In each Markov chain we generated A^]vic= 100,000 subsets. The subsets in the 
Mai^kov chain are coiTclated, producing A'mc/2t independent measurements for an integrated au- 
tocoiTclation time T. The statistical errors are determined taking the autocorrelations into account. 

The results of the subset method are compared with those computed with standard reweighting 
methods. For the latter we generate A'mc x N^ random matrices, such that the total number of 
matrices is the same as in the subset method. 
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Figure 2: Top row: chiral condensate E versus chemical potential jU^ for the subset method (blue bullets) and 
phase-quenched reweighting (red squares) for — 2,4,8. The full line corresponds to the exact analytical 
result [8]. Bottom row: relative statistical error e for the same data. The error for the reweighting method 
grows very rapidly and can only be trusted as long as the overlap problem is under control. 



Simulations were performed for = 2, . . . , 34 with 'Nf = 1 and m = 0.l /2N (the mass is small 
w.r.t. the magnitude of the smallest eigenvalue). In fig. 2 the results for the chiral condensate £ (top 
row) and its relative statistical error e (bottom row) are shown as a function of the chemical poten- 
tial. The statistical error of the phase-quenched reweighting grows exponentially with jU, until the 
method fails when the set of sampled matrices no longer overlaps with the relevant configurations. 
As the matrix size increases this failure occurs for smaller and smaller jj.^. This strongly contrasts 
with the subset method where the results are reliable up to much larger values of jJ.^ and agree with 
the analytical predictions of ref. [8]. Moreover, the error is independent of the chemical potential. 
We also investigated how the statistical error depends on the matrix size, see fig. 3. As expected, 
the work grows exponentially with A'^ for the reweighting method (right). However, for the subset 
method (left) the error is proportional to VN for a fixed number of sampled subsets, or A'^ for a 




Figure 3: Relative error e on the chiral condensate versus matrix size A^ for various values of the chemical 
potential. The results for the subset method are given on the left, those for phase-quenched reweighting on 
the right. Note that the later are given on a semi-log plot due to the exponential increase of the error. 
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fixed number of sampled matrices. Conversely, to achieve a constant error the number of sampled 
matrices should grow as N'^. 

A comparison of both methods makes clear how the sign problem is solved in the subset 
method: The cancellations needed to yield the exponentially small numbers in the partition function 
no longer happen through statistical sampling of the ensemble, but occur deterministically inside 
subsets of size of iff{N). Therefore, no exponential increase of the computing time with volume 
and chemical potential is expected in the subset method, as was confirmed by the numerical results. 

Note that in the reported study the fermionic weights (5.4) were directly computed by sum- 
ming over the complex determinants at chemical potential /i. However, we could equally well 
use eq. (5.5) to compute the exponentially small subset weights at /i from those at /i = 0, hence 
avoiding the need for any numerical cancellations. This strategy would be accompanied by some 
overhead because the determinants have to be computed both at = 0, to compute the subset 
weights, and at jU to compute the observable using eq. (6.1). 

8. Summary 

In this talk I have discussed the sign problem occurring in dynamical simulations of random 
matrices and presented a subset method which solves this problem for the Osborn model. The main 
feature of the method is the construction of subsets of matrices for which the fermionic weights, i.e. 
the sums of complex fermion determinants, ai^e real and positive. Importance sampling methods 
can be used to construct Markov chains of subsets and sample the random matrix ensemble. The 
numerical simulations confirmed that the subset method solves the sign problem for this model. 

The method relies on the ability to construct subsets with positive weights, and an important 
question arises as to what conditions need to be satisfied to apply this method to relevant physical 
systems. This matter will be explored further in future reseaixh. 
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